#******************************************************************************
#******************************************************************************

# standard

from dataclasses import dataclass, field

# from scipy.optimize import fsolve

# local

import numpy as np

from .twin import StandardisedTwinPipe

from .classes import StandardisedPipe

from ..common import formulas as core

from ..pipes import fic

from ..pipes import utils

#******************************************************************************
#******************************************************************************

HL_MODEL_DIRECT = 'direct'
HL_MODEL_SYM_ASYM = 'symasym'
    
TRGTPT_KRISCHER1936 = 'krischer1936'
TRGTPT_MULTIPOLE_ZERO_ORDER = 'wal1991mpzero'
TRGTPT_MULTIPOLE_FIRST_ORDER = 'wal1991mpfirst' 
TRGTPT_MULTIPOLE_SECOND_ORDER = 'wal1991mpsecond' 
TRGTPT_TWO_MODEL_APPROX = 'wal1991tmapprox'
    
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

# TODO: method to identify trench sizel limits for a given set of pipes

@dataclass
class SupplyReturnPipeTrenchWithIdenticalPipes:
    
    # provides pipe specific data that does not depend on the pipe length
    
    # constants
    pipes: list # list of StandardisedPipe objects for 1 line (supply or return)
    fluid_database: fic.FluidDatabase # should be only one
    ground_thermal_conductivity: float
    ground_air_heat_transfer_coefficient: float
    pipe_center_depth: list # can have the same size as the pipes list
    pipe_center_distance: list # not necessary in the case of twin pipes
    supply_temperature: float
    return_temperature: float
    max_specific_pressure_loss: float # in one pipe (supply or return)
    
    # time-varying parameters
    time_interval_duration: list # one element per time interval
    surroundings_temperature: list # one element per time interval
    
    # other parameters
    ground_thermal_resistance_model: str = TRGTPT_MULTIPOLE_ZERO_ORDER
    reference_heat_flow_rate: float = None
    reference_mass_flow_rate: list = field(init=False)
    reference_volumetric_flow_rate: list = field(init=False)
    reference_fluid_speed: list = field(init=False)
    nominal_fluid_speed: list = field(init=False) 
    nominal_volumetric_flow_rate: list = field(init=False) 
    nominal_mass_flow_rate: list = field(init=False) 
    nominal_heat_flow_rate: list = field(init=False) 
    heat_transfer_supply: list = field(init=False) # list of lists
    heat_transfer_return: list = field(init=False) # list of lists
    maximum_trench_length: list = field(init=False) # list of lists
    fluid_bulk: fic.Fluid = field(init=False)
               
    def __post_init__(self):
        
        # validate the inputs
        
        self.validate_inputs()
        
        # fluid properties
        
        self.fluid_bulk = fic.Fluid(
            phase='l',
            temperature=0.5*(self.supply_temperature+
                             self.return_temperature),
            pressure=1e5,
            db=self.fluid_database)
    
        # compute key performance parameters
        
        self.estimate_performance(
            self.ground_thermal_resistance_model)
        
        # validate outputs
        
    #**************************************************************************
    #**************************************************************************
    
    def validate_inputs(self):
        
        # validate the inputs
        
        # pipes can be a list of standardised pipe objects or a single one
        
        assert (
            type(self.pipes) == list or
            type(self.pipes) == StandardisedPipe or
            type(self.pipes) == StandardisedTwinPipe 
            )
        
        if type(self.pipes) == list:
            
            for pipe in self.pipes:
                
                #assert issubclass(type(pipe), StandardisedPipe)
                    
                assert (
                    type(pipe) == StandardisedPipe or
                    type(pipe) == StandardisedTwinPipe
                    )
                
        # fluid_database should be of type FluidDatabase
        
        assert type(self.fluid_database) == fic.FluidDatabase 
        
        # ground_air_heat_transfer_coefficient should be numeric
        
        assert (
            issubclass(
                type(self.ground_air_heat_transfer_coefficient), 
                float) 
            or
            issubclass(
                type(self.ground_air_heat_transfer_coefficient),
                int) 
            )
        
        # ground_thermal_conductivity should be numeric
            
        assert (
            issubclass(
                type(self.ground_thermal_conductivity), 
                float) 
            or
            issubclass(
                type(self.ground_thermal_conductivity),
                int) 
            )
        
        #**********************************************************************
        
        # pipe_center_depth can be numeric or a list
        
        assert (
            issubclass(
                type(self.pipe_center_depth), 
                float) 
            or
            issubclass(
                type(self.pipe_center_depth),
                int) 
            or
            type(self.pipe_center_depth) == list
            )
        
        # if it is a list, it needs to have the same size as self.pipes
        
        if type(self.pipe_center_depth) == list:
        
            assert len(self.pipes) == len(self.pipe_center_depth)
        
        #**********************************************************************
        
        # pipe_center_distance can be numeric or a list
        # note: it should be None for twin pipes since it is a pipe parameter
        
        assert (
            # single pipes: numeric or integer
            ((type(self.pipes) == StandardisedPipe 
              or 
              (type(self.pipes) == list and 
               (type(pipe) == StandardisedPipe for pipe in self.pipes))
              )
             and
             (issubclass(type(self.pipe_center_distance), float) or
              issubclass(type(self.pipe_center_distance), int) or
              type(self.pipe_center_distance) == list)
            )
            or
            # twin pipes: None
            ((type(self.pipes) == StandardisedTwinPipe or 
              (type(self.pipes) == list and
               (type(pipe) == StandardisedTwinPipe for pipe in self.pipes)))
             and
             self.pipe_center_distance == None)
            )
        
        # if it is a list, it needs to have the same size as self.pipes
        
        if type(self.pipe_center_distance) == list:
        
            assert len(self.pipes) == len(self.pipe_center_distance)
        
        #**********************************************************************
        
        # supply_temperature should be numeric
        
        assert (
            issubclass(
                type(self.supply_temperature), 
                float) 
            or
            issubclass(
                type(self.supply_temperature),
                int) 
            )
        
        # return_temperature should be numeric
        
        assert (
            issubclass(
                type(self.return_temperature), 
                float) 
            or
            issubclass(
                type(self.return_temperature),
                int) 
            )
        
        # max_specific_pressure_loss should be numeric
        
        assert (
            issubclass(
                type(self.max_specific_pressure_loss), 
                float) 
            or
            issubclass(
                type(self.max_specific_pressure_loss),
                int) 
            )
    
    #**************************************************************************
    #**************************************************************************
       
    def estimate_performance(
            self,
            ground_resistance_model=TRGTPT_MULTIPOLE_ZERO_ORDER):
        
        # initialise variables
        
        self.nominal_fluid_speed = []
        self.nominal_volumetric_flow_rate = []
        self.nominal_mass_flow_rate = []
        self.nominal_heat_flow_rate = []
        self.heat_transfer_supply = []
        self.heat_transfer_return = []
        self.maximum_trench_length = []
        
        self.reference_mass_flow_rate = []
        self.reference_volumetric_flow_rate = []
        self.reference_fluid_speed = []
        
        # run the calculations
        
        for i, pipe in enumerate(self.pipes):
            
            # self.pipe_center_depth can be a list or a number
            
            if type(self.pipe_center_depth) == list:
                
                pipe_center_depth = self.pipe_center_depth[i]
                
            else:
                
                pipe_center_depth = self.pipe_center_depth
                
            # self.pipe_center_distance can be a list, a number or None
            
            if isinstance(pipe, StandardisedTwinPipe):
                
                # the input should be None, the value should be extracted from
                # the pipe object
                
                pipe_center_distance = pipe.pipe_center_distance
                
            elif type(self.pipe_center_distance) == list:
                
                # list
                
                pipe_center_distance = self.pipe_center_distance[i]
                
            else:
                
                # numeric
                
                pipe_center_distance = self.pipe_center_distance
            
            # estimate performance for this pipe
            
            (re,
             f,
             nominal_fluid_speed,
             nominal_volumetric_flow_rate,
             nominal_mass_flow_rate,
             nominal_heat_flow_rate,
             list_q_s,
             list_q_r,
             list_ztel,
             reference_mass_flow_rate,
             reference_volumetric_flow_rate,
             reference_fluid_speed
             ) = self.estimate_pipe_performance(
                 pipe, 
                 pipe_center_depth, 
                 pipe_center_distance,
                 ground_resistance_model=ground_resistance_model)
                 
             # store results
             
            self.nominal_fluid_speed.append(
                nominal_fluid_speed)
            self.nominal_volumetric_flow_rate.append(
                nominal_volumetric_flow_rate)
            self.nominal_mass_flow_rate.append(
                nominal_mass_flow_rate)
            self.nominal_heat_flow_rate.append(
                nominal_heat_flow_rate)
            self.heat_transfer_supply.append(
                list_q_s)
            self.heat_transfer_return.append(
                list_q_r)
            self.maximum_trench_length.append(
                list_ztel)

            self.reference_mass_flow_rate.append(
                reference_mass_flow_rate)
            self.reference_volumetric_flow_rate.append(
                reference_volumetric_flow_rate)
            self.reference_fluid_speed.append(
                reference_fluid_speed)
        
    #**************************************************************************
    #**************************************************************************
        
    def estimate_pipe_performance(
            self, 
            pipe_obj: StandardisedPipe,
            pipe_center_depth: float,
            pipe_center_distance: float,
            ground_resistance_model=TRGTPT_MULTIPOLE_ZERO_ORDER) -> tuple:
        
        # calculates the pipe

        #**********************************************************************
        #**********************************************************************
        
        # What needs to be done?
        # 1) the nominal fluid speed needs to be determined
        # >> one value per pipe specification
        # 2) the design capacity needs to be determined
        # >> one value per pipe and network specification
        # 3) the specific heat losses need to be determined
        # >> one value per pipe, network and time interval
        
        #**********************************************************************
        #**********************************************************************
        
        # # fluid properties at bulk temperature
        
        # self.fluid_bulk = fic.Fluid(
        #     phase='l',
        #     temperature=0.5*(self.supply_temperature+
        #                      self.return_temperature),
        #     pressure=1e5,
        #     db=self.fluid_database)
        
        # self.nominal_fluid_speed = []
        # self.nominal_volumetric_flow_rate = []
        # self.nominal_mass_flow_rate = []
        # self.nominal_heat_flow_rate = []
        # self.heat_transfer_supply = []
        # self.heat_transfer_return = []
        
        # self.reference_mass_flow_rate = []
        # self.reference_volumetric_flow_rate = []
        # self.reference_fluid_speed = []
    
        #******************************************************************
        #******************************************************************
        
        # 1) the nominal fluid speed needs to be determined
    
        #******************************************************************
        #******************************************************************
        
        # for each pipe
        
        # generate initial solution
        
        if pipe_obj.e_eff == 0:
        
            # smooth pipes: the dff is a function of the reynolds number
             
            (x0_reynolds_number, 
             x0_friction_factor) = fic.DarcyFrictionFactorExplicit(
                  pipe_obj,
                  model=fic.DFF_PETHUKOV,
                  reynolds_number=None)
         
        else:
            
            # rough pipes: the dff is a function of the pipe roughness
              
            (x0_reynolds_number, 
             x0_friction_factor) = fic.DarcyFrictionFactorExplicit(
                  pipe_obj,
                  model=fic.DFF_NIKURADSE,
                  #reynolds_number=fic.DFF_RE_LIMITS[fic.DFF_NIKURADSE][1],
                  reynolds_number=None)
        
        # (x0_reynolds_number, 
        #  x0_friction_factor) = (
        #      utils.DarcyFrictionFactorReynoldsNumberForPipe(pipe_obj)
        #      )
                                                    
        x0_fluid_speed = core.MeanFluidSpeedFromReynoldsNumber(
            pipe_obj.d_int,
            x0_reynolds_number,
            self.fluid_bulk.kinematic_viscosity)
        
        #******************************************************************
        #******************************************************************
        
        # find the fluid speed
        
        for model in fic.DFF_MODELS:
            
            # try solving
        
            (v,
             re,
             f) = utils.find_specific_pressure_loss_matching_flow(
                pipe_obj,
                self.fluid_bulk,
                self.max_specific_pressure_loss, # per pipe
                dff_model=model,
                x0_fluid_speed=x0_fluid_speed,
                x0_friction_factor=x0_friction_factor,
                x0_reynolds_number=x0_reynolds_number)
                 
            if (v > 0 and re > 0 and f > 0):
                
                break
             
            #**************************************************************
            #**************************************************************
            
        nominal_fluid_speed = v
        
        #******************************************************************
        #******************************************************************
        
        # 2) the nominal capacity needs to be determined

        #******************************************************************
        #******************************************************************
        
        # volumetric flow rate
        
        nominal_volumetric_flow_rate = (
            core.VolumetricFlowRateFromMeanFluidSpeed(
                nominal_fluid_speed,
                pipe_obj.InternalCrossSectionalArea())
            )
        
        # mass flow rate
        
        nominal_mass_flow_rate = (
            core.MassFlowRateFromVolumetricFlowRate(
                nominal_volumetric_flow_rate, 
                self.fluid_bulk.mass_density)
            )
        
        # heat flow rate
        
        nominal_heat_flow_rate = (
            abs(core.HeatTransferRateMassFlow(
                self.supply_temperature,
                self.return_temperature,
                nominal_mass_flow_rate,
                self.fluid_bulk.specific_heat)
                )
            )
          
        #******************************************************************
        #******************************************************************
        
        # 3) the specific heat losses need to be determined
        
        #******************************************************************
        #******************************************************************
        
        # the reference fluid speed needs to be determined
        
        if self.reference_heat_flow_rate == None:
        
            # if the reference heat flow rate is None, use the nominal one
            
            reference_mass_flow_rate = (
                core.MassFlowRateFromHeatFlowRate(
                    nominal_heat_flow_rate,
                    abs(self.supply_temperature-self.return_temperature),
                    self.fluid_bulk.specific_heat)
                )
        
        elif self.reference_heat_flow_rate != 0: 
            
            # if not, and not zero either, use the reference heat flow rate
        
            reference_mass_flow_rate = (
                core.MassFlowRateFromHeatFlowRate(
                    self.reference_heat_flow_rate,
                    abs(self.supply_temperature-self.return_temperature),
                    self.fluid_bulk.specific_heat)
                )
            
        else: # reference heat flow rate is zero
        
            # and the reference mass flow rate is zero
        
            reference_mass_flow_rate = 0
            
        # internal heat transfer coefficient
        
        if reference_mass_flow_rate == 0:
            
            # if the mass flow rate is zero, the heat transfer coef. is inf
            
            reference_volumetric_flow_rate = 0
            
            reference_fluid_speed = 0
            
            internal_heat_transfer_coefficient = np.inf

        else: # if not zero
        
            # reference (not nominal) volumetric flow rate
            
            reference_volumetric_flow_rate = (
                core.VolumetricFlowRateFromMassFlowRate(
                    reference_mass_flow_rate,
                    self.fluid_bulk.mass_density)
                )
            
            # reference (not nominal) mean fluid speed
            
            reference_fluid_speed = (
                core.MeanFluidSpeedViaVolumetricFlowRate(
                    reference_volumetric_flow_rate,
                    pipe_obj.InternalCrossSectionalArea())
                )
                    
            # forced internal convection heat transfer coefficient
    
            internal_heat_transfer_coefficient = (
                HeatTransferCoefficientInternalForcedConvection(
                    pipe_obj,
                    self.fluid_bulk,
                    reference_fluid_speed)
                )
            
        #******************************************************************
        
        list_q_s = []
        
        list_q_r = []
        
        list_ztel = []
        
        for t in range(len(self.time_interval_duration)):
            
            # compute heat losses
            
            if isinstance(pipe_obj, StandardisedTwinPipe):
                    
                # twin pipes 
                
                (rsym,
                 rasym) = ThermalResistanceBuriedTwinPipesWallenten(
                    pipe_obj,
                    self.ground_thermal_conductivity,
                    pipe_center_depth,
                    pipe_center_distance,
                    self.ground_air_heat_transfer_coefficient,
                    model=ground_resistance_model
                    )
                     
                (q_s,
                 q_r) = HeatLossesSupplyReturnPipesSymmetricalAntiSymmetrical(
                    self.supply_temperature,
                    self.return_temperature,
                    self.surroundings_temperature[t],
                    rsym,
                    rasym)
            
            elif isinstance(pipe_obj, StandardisedPipe):
                
                # single pipes
                
                (u1,
                 u2,
                 _,
                 _,
                 _) = HeatTransferCoefficientsTwoBuriedSingles(
                      pipe_obj,
                      pipe_center_depth,
                      pipe_center_distance,
                      self.ground_thermal_conductivity,
                      self.ground_air_heat_transfer_coefficient,
                      internal_heat_transfer_coefficient,
                      model=ground_resistance_model
                     )
                     
                q_s, q_r = HeatLossesBuriedSupplyReturnPipesDirect(
                    self.supply_temperature,
                    self.return_temperature,
                    self.surroundings_temperature[t],
                    u1, 
                    u2)
            
            # else:
                
            #     raise NotImplementedError
                 
            # heat losses
            
            list_q_s.append(q_s)
             
            list_q_r.append(q_r)
            
            # zero efficiency trench length
            
            ztel = zero_transmission_efficiency_trench_length(
                design_trench_heat_transfer_rate=nominal_heat_flow_rate, 
                specific_trench_gains_or_losses=q_s+q_r,
                heating_network_rates_are_positive=( 
                    True if self.supply_temperature >= self.return_temperature
                    else False), # rates are always positive (absolute value)
                fluid_heat_losses_are_positive=True)
                        
            list_ztel.append(ztel)
            
        #**********************************************************************
        #**********************************************************************
        
        return (
            re,
            f,
            nominal_fluid_speed,
            nominal_volumetric_flow_rate,
            nominal_mass_flow_rate,
            nominal_heat_flow_rate,
            list_q_s,
            list_q_r,
            list_ztel,
            reference_mass_flow_rate,
            reference_volumetric_flow_rate,
            reference_fluid_speed
            )
    
    #**************************************************************************
    #**************************************************************************
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
    
def zero_transmission_efficiency_trench_length(
        design_trench_heat_transfer_rate: float,
        specific_trench_gains_or_losses: float,
        heating_network_rates_are_positive: bool = True,
        fluid_heat_losses_are_positive: bool = True):
    """
    Determines the trench length that effectively prevents heat transport.
    
    The calculation is dependent on the combined (supply and return) heat loss-
    es and on the design heat transfer rate. These two parameters can be provi-
    ded using different sign conventions. This calculation assumes that the 
    the specific heat losses or gains are constant along the pipe length.

    Parameters
    ----------
    design_trench_heat_transfer_rate : float
        The design heat transfer rate for the trench. The sign convention is
        determined by the parameter heating_network_rates_are_positive.
    specific_trench_gains_or_losses : float
        The specific (per trench meter) combined heat losses or gains for the 
        trench under consideration. The sign convention is determined by the
        parameter fluid_heat_losses_are_positive.
    heating_network_rates_are_positive : bool, optional
        If True, positive heat transfer rates imply a heating network. If not,
        they imply a cooling network. The default is True.
    fluid_heat_losses_are_positive : bool, optional
        If True, positive specific_trench_gains_or_losses values indicate heat
        losses. If False, they indicate heat gains. The default is True.

    Returns
    -------
    TYPE
        The trench length that prevents heat transport.

    """
    
    #**************************************************************************
    #**************************************************************************
    
    # two cases have valid finite lengths:
    # 1) heating network with heat losses
    # 2) cooling network with heat gains
    # other cases have infinite lengths:
    # 1) heating network with heat gains
    # 2) cooling network with heat losses
    
    if (heating_network_rates_are_positive and 
        fluid_heat_losses_are_positive):
        
        # polarities: 
        # i) design rates are for heating if positive, for cooling if negative
        # ii) specific losses are positive, specific gains are negative
        
        return (
            design_trench_heat_transfer_rate/specific_trench_gains_or_losses
            if ((design_trench_heat_transfer_rate >= 0 and 
                 specific_trench_gains_or_losses > 0) or
                (design_trench_heat_transfer_rate <= 0 and
                 specific_trench_gains_or_losses < 0))
            else np.inf
            )
    
    elif (heating_network_rates_are_positive and 
          not fluid_heat_losses_are_positive):
          
          # polarities: 
          # i) design rates are for heating if positive, for cooling if negative
          # ii) specific losses are negative, specific gains are positive
      
          return (
              -design_trench_heat_transfer_rate/specific_trench_gains_or_losses
              if ((design_trench_heat_transfer_rate >= 0 and 
                   specific_trench_gains_or_losses < 0) or
                  (design_trench_heat_transfer_rate <= 0 and
                   specific_trench_gains_or_losses > 0))
              else np.inf
              )
    
    elif (not heating_network_rates_are_positive and 
          fluid_heat_losses_are_positive):
          
          # polarities: 
          # i) design rates are for cooling if positive, for heating if negative
          # ii) specific losses are positive, specific gains are negative
      
          return (
              -design_trench_heat_transfer_rate/specific_trench_gains_or_losses
              if ((design_trench_heat_transfer_rate <= 0 and 
                   specific_trench_gains_or_losses > 0) or
                  (design_trench_heat_transfer_rate >= 0 and
                   specific_trench_gains_or_losses < 0))
              else np.inf
              )
    
    # elif (not heating_network_rates_are_positive and 
    #       not fluid_heat_losses_are_positive):
    
    else: 
    
        # polarities: 
        # i) design rates are for cooling if positive, for heating if negative
        # ii) specific losses are negative, specific gains are positive
        
        return (
            design_trench_heat_transfer_rate/specific_trench_gains_or_losses
            if ((design_trench_heat_transfer_rate <= 0 and 
                 specific_trench_gains_or_losses < 0) or
                (design_trench_heat_transfer_rate >= 0 and
                 specific_trench_gains_or_losses > 0))
            else np.inf
            )

    #**************************************************************************
    #**************************************************************************
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def HeatTransferCoefficientsTwoBuriedSingles(
     pipe,
     pipe_center_depth,
     pipe_center_distance,
     ground_thermal_conductivity,
     ground_air_heat_transfer_coefficient,
     internal_heat_transfer_coefficient,
     model=TRGTPT_MULTIPOLE_FIRST_ORDER):
     
    # compute the heat losses using the Bohm (2000) model
    
    #**************************************************************************
    
    # parameters
    
    #beta = np.log(pipe.d_ins/pipe.d_ext)*ground_thermal_conductivity/pipe.k_ins
    
    # modified pipe depth to account for ground-air heat transfer
    
    pipe_center_depth_modified = ( 
        pipe_center_depth+(ground_thermal_conductivity/
                           ground_air_heat_transfer_coefficient)
        )
    
    # resistances
    
    # resistance of the ground (ground and pipe casing surface interfaces too)
        
    r_g = ThermalResistanceGroundTwoPipeTrench(
        pipe,
        ground_thermal_conductivity,
        ground_air_heat_transfer_coefficient,
        pipe_center_depth,
        pipe_center_distance,
        model=model)
    
    # resistance across the insulation and casing
    
    r_i = fic.InsulatedPipeSpecificResistance(
        h_internal=internal_heat_transfer_coefficient,
        h_external=np.inf, 
        extra_specific_resistance=0, 
        pipe=pipe)
    
    # resistance relative to the interaction between the pipes
    
    r_m = (
        np.log(1+(2*pipe_center_depth_modified/pipe_center_distance)**2)/
        (4*np.pi*ground_thermal_conductivity)
        )
    
    # u1 and u2 coefficients
    
    u1 = (r_g+r_i)/((r_g+r_i)**2-r_m**2)
    
    u2 = r_m/((r_g+r_i)**2-r_m**2)
    
    # return
    
    return u1, u2, r_g, r_i, r_m
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def HeatLossesBuriedSupplyReturnPipesDirect(temperature_supply,
                                           temperature_return,
                                           temperature_ground,
                                           u1,
                                           u2):
    
    # rsym = u1-u2
    # rasym = u1+u2
    
    q_s = (
        (u1-u2)*(temperature_supply-temperature_ground)
        +
        u2*(temperature_supply-temperature_return)
        )
        
    q_r = (
        (u1-u2)*(temperature_return-temperature_ground)
        -
        u2*(temperature_supply-temperature_return)
        )
        
    # if model == 0:

    #     # Simpler alternative:
                                             
    #     q_s = (
    #         (u1-u2)*(temperature_supply-temperature_ground)
    #         +
    #         u2*(temperature_supply-temperature_return)
    #         )
            
    #     q_r = (
    #         (u1-u2)*(temperature_return-temperature_ground)
    #         -
    #         u2*(temperature_supply-temperature_return)
    #         )
    
    # else:
    
    #     # Wallenten, 1991
    
    #     q_s = (
    #         (u1-u2)*
    #         ((temperature_supply+temperature_return)/2-temperature_ground)
    #         +
    #         (u1+u2)*(temperature_supply-temperature_return)/2
    #         )
            
    #     q_r = (
    #         (u1-u2)*
    #         ((temperature_supply+temperature_return)/2-temperature_ground)
    #         -
    #         (u1+u2)*(temperature_supply-temperature_return)/2
    #         )
    
    return q_s, q_r
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def ThermalResistanceGroundTwoPipeTrench(pipe,
                                         ground_thermal_conductivity,
                                         ground_air_heat_transfer_coefficient,
                                         pipe_center_depth,
                                         pipe_center_distance,
                                         model=TRGTPT_MULTIPOLE_FIRST_ORDER):
        
    beta = np.log(pipe.d_ins/pipe.d_ext)*ground_thermal_conductivity/pipe.k_ins
    
    # modified pipe depth to account for ground-air heat transfer
    
    pipe_center_depth_modified = ( 
        pipe_center_depth+(ground_thermal_conductivity/
                           ground_air_heat_transfer_coefficient)
        )
    
    if model == TRGTPT_KRISCHER1936:
        
        return (
            np.log(2*pipe_center_depth_modified/pipe.d_cas+
                    np.sqrt((2*pipe_center_depth_modified/pipe.d_cas)**2-1))/
            (2*np.pi*ground_thermal_conductivity)
            )
        
    elif model == TRGTPT_MULTIPOLE_ZERO_ORDER:
        
        return (
            np.log(4*pipe_center_depth_modified/pipe.d_cas)/
            (2*np.pi*ground_thermal_conductivity)
            )
        
    elif model == TRGTPT_MULTIPOLE_FIRST_ORDER:
    
        et_numerator = (
            (pipe.d_cas/(2*pipe_center_distance))**2
            +
            (pipe.d_cas/(4*pipe_center_depth_modified))**2
            +
            (pipe.d_cas**2)/(4*(4*pipe_center_depth_modified**2+
                                pipe_center_distance**2))
            )
        
        et_denominator = (
            (1+beta)/(1-beta)
            -
            (pipe.d_cas/(2*pipe_center_distance))**2
            )
            
        return (
            (np.log(4*pipe_center_depth_modified/pipe.d_cas)
              -
              et_numerator/et_denominator)/
            (2*np.pi*ground_thermal_conductivity)
            )
        
    else:
        
        raise NotImplementedError

#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def HeatTransferCoefficientInternalForcedConvection(
        pipe,
        fluid,
        fluid_speed):
    
    # reynolds number
    
    reynolds_number = core.ReynoldsNumber(pipe.d_int,
                                         fluid_speed,
                                         fluid.kinematic_viscosity)
    
    # flow regime
    
    flow_regime = fic.FlowRegime(reynolds_number)
    
    # friction factor
    
    friction_factor = fic.DarcyFrictionFactor(reynolds_number,
                                              pipe,
                                              flow_regime)
    
    # nusselt number
    
    nusselt_number = fic.NusseltNumber(pipe,
                                       friction_factor,
                                       reynolds_number,
                                       fluid.prandtl_number,
                                       flow_regime)
    
    # heat transfer coefficient
    
    return core.ConvectionHeatTransferCoefficient(pipe.d_int,
                                                  nusselt_number,
                                                  fluid.thermal_conductivity)
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def ThermalResistanceBuriedSinglePipeWallenten(
        pipe,
        soil_k,
        pipe_depth,
        soil_h,
        model=TRGTPT_MULTIPOLE_SECOND_ORDER):
    
    beta = (soil_k/pipe.k_ins)*np.log(pipe.d_cas/pipe.d_ext)
    
    adjusted_pipe_depth = pipe_depth+soil_k/soil_h
    
    if model == TRGTPT_MULTIPOLE_ZERO_ORDER:
        
        return (
            (np.log(4*adjusted_pipe_depth/pipe.d_cas)
             +
             beta)/
            (2*np.pi*soil_k)
            )
    
    elif model == TRGTPT_MULTIPOLE_FIRST_ORDER:
        
        return (
            (np.log(4*adjusted_pipe_depth/pipe.d_cas)
             +
             beta
             +
             1/(1-((1+beta)/(1-beta))*(4*adjusted_pipe_depth/pipe.d_cas)**2))/
            (2*np.pi*soil_k)
            )
    
    elif model == TRGTPT_MULTIPOLE_SECOND_ORDER:
        
        return (
            (np.log(4*adjusted_pipe_depth/pipe.d_cas)
             +
             beta
             +
             1/(1-((1+beta)/(1-beta))*(4*adjusted_pipe_depth/pipe.d_cas)**2)
             +
             (1
              +
              ((1+beta)*(1-2*beta)*(pipe.d_cas/(4*adjusted_pipe_depth))**2)/
              (2*(1-beta)*(1+2*beta))
              -
              (3*(1-2*beta)*(pipe.d_cas/(4*adjusted_pipe_depth))**4)/
              (2*(1+2*beta))
              )
             /
             (1
              -
              ((4*adjusted_pipe_depth/pipe.d_cas)**2-
               3*(4*adjusted_pipe_depth/pipe.d_cas)**2)*(1+beta)/(1-beta)
              -
              ((pipe.d_cas/(4*adjusted_pipe_depth))**4)*(1-2*beta)/(1+2*beta)
              )
             )/
            (2*np.pi*soil_k)
            )
    
    else:
        
        raise NotImplementedError 
            
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def ThermalResistanceTwoBuriedSinglePipesWallenten(
        pipe,
        soil_k,
        pipe_depth,
        pipe_distance,
        soil_h,
        model=TRGTPT_MULTIPOLE_FIRST_ORDER):
    
    beta = (soil_k/pipe.k_ins)*np.log(pipe.d_cas/pipe.d_ext)
    
    adjusted_pipe_depth = pipe_depth+soil_k/soil_h
    
    if model == TRGTPT_MULTIPOLE_ZERO_ORDER:
        
        r_s = (
            (np.log(4*adjusted_pipe_depth/pipe.d_cas)
             +
             beta
             +
             np.log(np.sqrt(1+(2*adjusted_pipe_depth/pipe_distance)**2)))/
            (2*np.pi*soil_k)
            )
                
        r_a = (
            (np.log(4*adjusted_pipe_depth/pipe.d_cas)
             +
             beta
             -
             np.log(np.sqrt(1+(2*adjusted_pipe_depth/pipe_distance)**2)))/
            (2*np.pi*soil_k)
            )
    
    elif model == TRGTPT_MULTIPOLE_FIRST_ORDER:
        
        r_s = (
            (np.log(4*adjusted_pipe_depth/pipe.d_cas)
             +
             beta
             +
             np.log(np.sqrt(1+(2*adjusted_pipe_depth/pipe_distance)**2))
             -
             ((pipe.d_cas/(2*pipe_distance))**2
              +
              (pipe.d_cas/(4*adjusted_pipe_depth))**2
              +
              (1/16)*(pipe.d_cas**2)/
              ((pipe_distance/2)**2+(adjusted_pipe_depth)**2))/
             ((1+beta)/(1-beta)+(pipe.d_cas/(2*pipe_distance))**2))/
            (2*np.pi*soil_k)
            )
                
        r_a = (
            (np.log(4*adjusted_pipe_depth/pipe.d_cas)
             +
             beta
             -
             np.log(np.sqrt(1+(2*adjusted_pipe_depth/pipe_distance)**2))
             -
             ((pipe.d_cas/(2*pipe_distance))**2
              +
              (pipe.d_cas/(4*adjusted_pipe_depth))**2
              -
              (3/16)*(pipe.d_cas**2)/
              ((pipe_distance/2)**2+(adjusted_pipe_depth)**2))/
             ((1+beta)/(1-beta)-(pipe.d_cas/(2*pipe_distance))**2))/
            (2*np.pi*soil_k)
            )
    
    else:
        
        raise NotImplementedError
    
    return r_s, r_a
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def ThermalResistanceBuriedTwinPipesWallenten(
        pipe,
        soil_k,
        pipe_depth,       # H
        pipe_distance,    # 2*D
        soil_h,
        model=TRGTPT_MULTIPOLE_ZERO_ORDER):
    
    # TODO: redo using the original work's nomenclature
    
    #**************************************************************************

    sigma = (pipe.k_ins-soil_k)/(pipe.k_ins+soil_k)
    
    adjusted_pipe_depth = pipe_depth+soil_k/soil_h
    
    if model == TRGTPT_MULTIPOLE_ZERO_ORDER:
    
        # symmetrical
        
        r_s = (
            ((2*pipe.k_ins/soil_k)*np.log(4*adjusted_pipe_depth/pipe.d_cas)
             +
             np.log(((pipe.d_cas/2)**2)/(pipe_distance*pipe.d_ext/2))
             +
             sigma*
             np.log(
                 ((pipe.d_cas/2)**4)/
                 ((pipe.d_cas/2)**4-(pipe_distance/2)**4)
                 )
             )/
            (2*np.pi*pipe.k_ins)
            )
                
        # asymmetrical
        
        r_a = (
            (np.log(2*pipe_distance/pipe.d_ext)
             +
             sigma*
             np.log(
                 ((pipe.d_cas/2)**2+(pipe_distance/2)**2)/
                 ((pipe.d_cas/2)**2-(pipe_distance/2)**2))
             )/
            (2*np.pi*pipe.k_ins)
            )
                
    elif model == TRGTPT_MULTIPOLE_FIRST_ORDER:
        
        #**********************************************************************
        
        # numerator for the rs's extra term
    
        rs_et_num = (
            (pipe.d_ext/(2*pipe_distance))
            -
            (sigma*pipe.d_ext*(pipe_distance/2)**3)/
            ((pipe.d_cas/2)**4-(pipe_distance/2)**4)
            )**2
        
        # denominator for the rs's extra term
        
        rs_et_den = (
            1+(pipe.d_ext/(2*pipe_distance))
            +
            sigma*((pipe_distance*(pipe.d_ext/2)*(pipe.d_cas/2)**2
                    /((pipe.d_cas/2)**4-(pipe_distance/2)**4))**2)
            )
                
        r_s = (
            ((2*pipe.k_ins/soil_k)*np.log(4*adjusted_pipe_depth/pipe.d_cas)
             +
             np.log(((pipe.d_cas/2)**2)/(pipe_distance*pipe.d_ext/2))
             +
             np.log(
                 ((pipe.d_cas/2)**4)/
                 ((pipe.d_cas/2)**4-(pipe_distance/2)**4)
                 )*sigma
             -
             rs_et_num/rs_et_den)/
            (2*np.pi*pipe.k_ins)
            )
        
        #**********************************************************************
            
        gamma = (
            2*(1-sigma**2)/(1-sigma*(pipe.d_cas/(4*adjusted_pipe_depth))**2)
            )
        # gamma = (
        #     -2*(1+sigma)/(1-sigma*(pipe.d_cas/(4*adjusted_pipe_depth))**2)
        #     )
        
        # TODO: derive the equations from scratch
        
        ra_et_num = (
            pipe.d_ext/(2*pipe_distance)
            -
            gamma*(pipe_distance/2)*(pipe.d_ext/2)/(4*adjusted_pipe_depth**2)
            -
            pipe_distance*sigma*(pipe.d_ext/2)*((pipe.d_cas/2)**2)/
            ((pipe.d_cas/2)**4-(pipe_distance/2)**4)
            )**2
            
        ra_et_den = (
            1-(pipe.d_ext/(2*pipe_distance))**2
            -
            gamma*(pipe.d_ext/2)/(2*adjusted_pipe_depth)
            +
            2*sigma*((pipe.d_ext/2)**2)*((pipe.d_cas/2)**2)*
            ((pipe.d_cas/2)**4+(pipe_distance/2)**4)/
            (((pipe.d_cas/2)**4-(pipe_distance/2)**4)**2)
            )
        
        r_a = (
            (np.log(2*pipe_distance/pipe.d_ext)
             +
             np.log(
                 ((pipe.d_cas/2)**2+(pipe_distance/2)**2)/
                 ((pipe.d_cas/2)**2-(pipe_distance/2)**2))*
             sigma
             -
             ra_et_num/ra_et_den
             -
             gamma*(pipe_distance/(4*adjusted_pipe_depth))**2)/
            (2*np.pi*pipe.k_ins)
            )
        
        #**********************************************************************
            
    elif model == TRGTPT_TWO_MODEL_APPROX:
        
        R1 = np.arccosh(0.5*(
            pipe.d_ext/pipe.d_cas+
            pipe.d_cas/pipe.d_ext-
            (pipe.d_cas/pipe.d_ext)*
            (pipe_distance/pipe.d_cas)**2)
            )
        
        R2 = 4*np.arccosh(2*(pipe_distance/pipe.d_ext)**2-1)
        
        R3 = (2*pipe.k_ins/soil_k)*np.log(
            2*adjusted_pipe_depth/pipe.d_cas+
            np.sqrt(
                (2*adjusted_pipe_depth/pipe.d_cas)**2
                -1)
            )
        
        r_s = (
            (R3+
             1/(1/R1
                -
                1/R2))/
            (2*np.pi*pipe.k_ins)
            )
        
        r_a = np.inf
        
    else:
        
        raise NotImplementedError

    #**************************************************************************
        
    # return
    
    return r_s, r_a
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def HeatLossesSupplyReturnPipesSymmetricalAntiSymmetrical(supply_temperature,
                                                          return_temperature,
                                                          ground_temperature,
                                                          rsym,
                                                          rasym):

    #**************************************************************************
    
    qs = (
        ((supply_temperature+return_temperature)/2-ground_temperature)/rsym
        +
        ((supply_temperature-return_temperature)/2)/rasym
        )
    
    qr = (
        ((supply_temperature+return_temperature)/2-ground_temperature)/rsym
        -
        ((supply_temperature-return_temperature)/2)/rasym
        )

    #**************************************************************************
    
    return qs, qr
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************

def HeatTransferCoefficientsFromSymmetricalAsymmetricalResistances(r_sym,
                                                                   r_asym):
    
    # u1122 = 0.5*(1/r_s+1/r_a)
    
    # u12 = 0.5*(1/r_a-1/r_s)
    
    return 0.5*(1/r_sym+1/r_asym), 0.5*(1/r_asym-1/r_sym)
        
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************
#******************************************************************************